Accurate Thermodynamics for Short-Ranged Truncations of Coulomb Interactions in 

Site-Site Molecular Models 
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Coulomb interactions are present in a wide variety of all-atom force fields. Spherical truncations 
of these interactions permit fast simulations but are problematic due to their incorrect thermo- 
dynamics. Herein we demonstrate that simple analytical corrections for the thermodynamics of 
uniform truncated systems are possible. In particular results for the SPC/E water model treated 
with spherically-truncated Coulomb interactions suggested by local molecular field theory [Proc. 
Nat. Acad. Sci. USA 105, 19136 (2008)] are presented. We extend results developed by Chan- 
dler [J. Chem. Phys. 65, 2925 (1976)] so that we may treat the thermodynamics of mixtures of 
flexible charged and uncharged molecules simulated with spherical truncations. We show that the 
energy and pressure of spherically-truncated bulk SPC/E water are easily corrected using exact 
second-moment-like conditions on long-ranged structure. Furthermore, applying the pressure cor- 
rection as an external pressure removes the density errors observed by other research groups in NPT 
simulations of spherically-truncated bulk species. 



I. INTRODUCTION 

Most classical intermolecular potential models assign 
effective point charges to intramolecular interaction sites 
to describe charge separation in polar molecules and the 
ability to form hydrogen bonds^ 2 -. Thus, even in purely 
neutral systems, charge-charge interactions remain im- 
portant and expensive components of molecular simula- 
tions, usually dealt with via Ewald summations or some 
other lattice-sum- like technique 3 . 

Recently there has been renewed interest in spherically 
truncating the 1/r interaction and neglecting the long- 
ranged components beyond a specified cutoff radius^^i 
This permits fast and efficient simulations that scale lin- 
early with system size. However, spherical truncation 
are problematic to implement for Coulomb interactions. 
While many groups have found that accurate local pair 
correlation functions in uniform systems may be obtained 
by a variety of spherical truncations of 1/r— i 5 -! 6 -^, two 
common and valid objections to spherical truncations re- 



1. they fail for structural and electrostatic properties 
in nonuniform systems, e.g., systems with point 
charges confined between wall s 9 ' 10 , and 

2. in uniform systems, the thermodynamics predicted 
by such truncations^ and even the bulk densities in 
NPT simulation s 7 ' 11 are known to be inaccurate. 

Recently we overcame the first objection, showing that 
local molecular field (LMF) theor y 12 ' 13 provides an accu- 
rate path to structural properties in both ionic and aque- 
ous nonuniform systems using a spherical truncation of 
1/r along with a restructured external potential Vr to 
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FIG. 1: Plot of total potential energy without (Uo/N in red 
crosses) and with the long-range correction (Uo/N + Ui/N 
in blue circles) for the full range of a studied, representing 
greater inclusion of nearby core interactions. The length a 
sets the scale for the smooth truncation of the Coulomb in- 
teractions as is explained further in Section [TTJ Error bars are 
smaller than the data points. The Ewald determined energy 
is indicated by a horizontal line. 



account for the net averaged effects of the long-ranged 
forces neglected in the spherical truncatio n 14 ' 15 ' 16 . Elec- 
trostatic properties are then also very accurately de- 
scribedi 6 -^. 

Generalizing previous work for purely ionic sys- 
tem s 18 ' 19 , we show here that the LMF framework also 
guides us to simple analytic corrections for the energy 
and pressure of a general uniform mixture of both polar 
and charged site-site molecules. As in thermodynamic 
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perturbation theory 2 ^, we can view a short-ranged trun- 
cation of 1 jr as the reference system for the fully interact- 
ing system. Our corrections are appropriate and accurate 
only for certain special well-chosen reference systems as 
discussed below, which we refer to as "mimic systems". 
The total energy and pressure of the full system is then 
given by the sums 

Utot = U + Ui 

Ptot=P + Pl- (1) 

Typical simulations using such a spherical truncation of 
the 1/r interaction yield only Uq and Pq, and our task is 
to determine the corrections U\ and P\ that would arise 
from an accurate treatment of the long-ranged interac- 
tions. 

In particular, as shown in Fig. Q]for bulk SPC/E wa- 
ter—, Uq alone over a wide range of truncation distances 
parametrized by the length a does not agree with the en- 
ergy as calculated using three-dimensional Ewald sums. 
In contrast, any a of 3.0 A or greater reproduces the 
short-ranged pair correlations predicted by Ewald sums 
quite wel l 16 ' 22 . As a increases, the energies are in better 
agreement with the Ewald calculated values, but notice- 
able differences remain even for large a. 

Related problems arose in recent NPT ensemble simu- 
lations of water— Using Wolf sums 5 -, researchers simu- 
lated systems quite similar to the truncated system dic- 
tated by LMF theory for a « 5.0 A (labeled "DFS 2 " in 
Ref. 0) and for a range of a in Ref. [111. Each used the 
NPT ensemble with a pressure of 1 atm. They found 
generally good structural agreement, but noted thermo- 
dynamic discrepancies, like an elevated energy and a de- 
pressed density as compared to Ewald simulations. 

Here we use LMF theory along with long-wavelength 
constraints on the behavior of charge correlation func- 
tions implied by exact expressions for the dielectric con- 
stant in neutral systems 2 ^ and the related Stillinger- 
Lovctt moment conditions 24 for ionic systems to derive 
analytic expressions for U\ and P\. Results incorporat- 
ing this correction for the energy of SPC/E water are also 
given in Fig. [T] and their high accuracy is evident. 



II. LOCAL MOLECULAR FIELD THEORY FOR 
SITE-SITE MOLECULAR SIMULATIONS 

Local molecular field theory provide a general theo- 
retical framework for assessing and correcting spherical 
truncations of Coulomb interactions. The derivation of 
LMF theory for systems with Coulomb interactions has 
been recently reviewed elsewhere^, and we will be brief 
in our discussion here. 

LMF theory divides the 1/r potential into short- and 
long-ranged parts characterized by the length a as 

1 , . , . erfc(r/cr) erf(r/a) 

-=v (r)+ Vl (r = K -L-L + V 1 ' . 2 

r r r 



This potential separation isolates strong short-ranged 
and rapidly- varying Coulomb interactions in t>o(r), while 
the remaining slowly- varying long-ranged forces are con- 
tained in v\(r). vi(r) is proportional to the electro- 
static potential arising from a smooth normalized Gaus- 
sian charge distribution with width a, and is defined by 
the convolution 

By construction, «i(r) is slowly- varying in r-space over 
the smoothing length a (see, e.g., Fig. 1 in Ref. Il3f). and 
contains only small wave vectors in reciprocal space, as 
can be seen from its Fourier transform 

4-TT 

v 1 (k) = ^eM~(ka) 2 /^}. (4) 

The short-ranged «o(r) = l/r—vx(r) is then the screened 
potential resulting from a point charge surrounded by a 
neutralizing Gaussian charge distribution whose width a 
also sets the scale for the smooth truncation of vo . At dis- 
tances much less than a the force from «o(r) approaches 
that from the full 1/r potential. 

Starting from the exact Yvon-Born-Green hierarchy 2 ^, 
and exploiting the slowly- varying nature of v\ (r) , LMF 
theory accounts for the averaged effects of the long- 
ranged component vi(r) in a mean-field sense by a 
rescaled, self-consistent, mean electrostatic potential 
Vr(t). The short-ranged vo(r), with a chosen large 
enough to capture relevant nearest-neighbor interactions 
like core repulsions and hydrogen bonding, is the spher- 
ical truncation used in LMF theory. A system in the 
presence of Vr(i") with 1/r replaced by the short-ranged 
VQ(r) is referred to as a mimic system, and densities as- 
sociated with such a system are indicated by PR(r). 

Often for uniform systems, Vr(i") has negligible ef- 
fect on short-ranged pair correlations 18,22 . For such uni- 
form systems, we may simulate simply using vq(t) with 
Vr(i") = and thus generate densities po( r )- We call this 
approximation to the full LMF theory the strong- coupling 
approximation and refer to the resulting truncated water 
model as Gaussian-truncated water. While not generally 
true, LMF theory in the strong-coupling approximation 
is related to other spherical truncations such as site-site 
reaction field 8 and Wolf truncations 5 . 

LMF theory seeks to obtain the properties of the 
full, uniform system from the simulation of the short- 
ranged system, whose total energy Uq for Gaussian- 
truncated SPC/E water includes all contributions from 
the Lennard- Jones interactions as well as the short- 
ranged components of the Coulomb interactions due to 
vo(r). Other authors have used the vo(r) truncation and 
proposed numerical corrections to the energy and pres- 
sure of ionic systems based on integral equation meth- 
ods^ 5 -, but the simple and accurate analytical corrections 
possible using moment conditions and our choice of v q (r) 
and v i (r) as described below in Sections IIVI and [V] have 
not been previously derived. 
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FIG. 2: Diagram of Gaussian-truncated SPC/E water. The 
traditional SPC/E model of water is shown in the center with 
three point charges and a Lennard- Jones core to represent the 
excluded volume. Gaussian-truncated water is constructed 
by replacing the three point charges with the corresponding 
short-ranged vo(r). 



III. SIMULATION DETAILS 

The main result of this paper is a general derivation of 
simple analytical corrections for the spherical truncation 
of Coulomb interactions in simple charged and uncharged 
site-site molecular models. In order to demonstrate the 
accuracy of these corrections, we also have carried out a 
series of simulations of a molecular water model at am- 
bient conditions. 

The water model we choose is SPC/E water—, shown 
in the center of Fig. [2] A Lennard- Jones core, depicted by 
the solid circle with diameter olj = 3.161 A accounts for 
the excluded volume of the molecules, and point charges 
are present at each of the atomic sites in order to rep- 
resent the charge separation along the OH bonds and to 
allow for hydrogen bonding between molecules. In order 
to simulate Gaussian-truncated water, we replace the 1/r 
interaction from each of these point charges by the short- 
ranged vo(r) as represented by the dashed circles drawn 
here to scale with diameter a = 4.5 A. 

We carried out a molecular dynamics simulation of a 
uniform system of 1728 SPC/E water molecules using a 
modified version of DLPOLy2.162£. The Berendsen ther- 
mostat^ with a relaxation time constant of 0.5 ps is used 
to maintain the temperature at 300 K, and, for the fi- 
nal set of data presented, a Berendsen barostat main- 
tains the pressure. All simulations use a timestep of l fs. 
For the spherical truncations, a ranges from 3.0 A to 
6.0 A, with the cutoff radius ranging from 9.5 A (the 
cutoff radius for the Lennard- Jones core ) to 13.5 A. As 
explored in Refs. [l6| and [2^, any a of 3.0 A or greater 
reproduces the short-ranged pair correlations predicted 
by Ewald sums quite well. For our benchmark, we com- 
pare to simulations using three-dimensional Ewald sums 
with a = 0.30 A -1 and fc max = 10. The systems were 



each equilibrated for a total of 500 ps, with 1.5 ns of data 
collection; error bars were based on 100 ps blocks of data. 

IV. ANALYTICAL ENERGY CORRECTION 
VIA MOMENT CONDITIONS 

We assume here that charged interactions arise only 
between charges on different molecules, as is the case for 
typical molecular liquid models like SPC/E water. These 
ideas can be extended to larger molecular species with in- 
tramolecular charge-charge interactions between further- 
neighbor sites, as briefly discussed in Appendix lAl 

The total Coulomb energy U q for the full system 
can then be exactly expressed in terms of a two-point 
intermolecular charge-density function p qq (with units 
charge 2 /volume 2 ) as 

2 J J |r-r'| 

= Y- J drp qq (r)v a (r) +Y J dv p qq (r) Vl {r), (5) 

where we have used the uniformity of the fluid and Eq. 
([2]) in the second equality. The composite function p qq (r) 
is a charge- weighted linear combination of all intermolec- 
ular, two-point site-site distribution function o 23 ' 28 , and a 
detailed expression is given in Appendix [A"l 

As in Ref. 0, we argue that the first term on the right 
in Eq. ([5]) can be accurately approximated by Uq, the en- 
ergy obtained directly from the Gaussian-truncated wa- 
ter simulation using vo(r) alone, because at short dis- 
tances where vq(t) is non- negligible, Po q (r) closely re- 
sembles the exact p qq [r). Thus we have 

^r*\J drp qq (r)v (r) + \j dr p qq '(r)«i(r) . (6) 

However, as noted in Ref. Il8l . a similar approximation 
for the second term will fail because «i(r) mainly con- 
tains small- wavevector components, exactly the range of 
fc-components where p\^ (r) will not accurately represent 
the p qq (r) of the full system. In fact, this integral will 
diverge if constraints due to neutrality in ionic systems 
are not obeyed. Similar considerations are true for the 
mixed molecular systems considered here. 

Thus, we again follow the more fruitful path of writ- 
ing the second term in fc-space and approximating the 
long wavelength behavior of the charge energy function 
based on exact relations. For a uniform system, we may 
use Parseval's relation and Eq. ([4| to reexpress Eq. ([6]) 
exactly as 

El ~ Tji + I_L_ f dk l^m e -fc 2 - 2 /4 (7) 

This choice in Eq. ([4]) of V\{r) in LMF theory allows 
us to make a highly useful approximation that high- 
lights its advantages over other possible potential sep- 
arations. The Gaussian from v\(k) damps out the large 
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fc-contributions to the second term, U q . Thus for suffi- 
ciently large a we can simply represent p qq (k) by its two 
smallest moments in fc-space, 



p qq (k) S3 p^ qq + p^ qq k 2 + 0(k A ). 



(8) 



We show in Appendix [Al that p qq (k) is very simply re- 
lated to the basic charge-charge linear response function 
X qq (k) that appears in Chandler's formula 23 for the di- 
electric constant of a neutral molecular mixture. We have 
generalized the derivation to include both neutral and 
charged molecular species in Appendix [X] and demon- 
strate the simpler expansion of x qq (k) in Appendix [B] 
This expression is a consequence of Stillinger-Lovett-like 
sum rules 20,2 - arising from the assumption that the po- 
tential induced by a test charge Q in a uniform molecular 
fluid approaches AnQ/ek 2 at small k to linear order in Q. 
This allows us to relate the moments of p qq to the dielec- 
tric constant e and other molecular properties. A system 
with mobile ions exhibits complete screening with e = oo. 

Here we simply state the final expansion of the two- 
point charge-density function up to second order in fc, 
as derived in Appendix [SJ We find for a general mix- 
ture of charged (C) and neutral (N) polarizable site-site 
molecules without intramolecular charge-charge interac- 
tions, 



P qq (k) = -J2pcq 2 c + k 



2 k B T 6-1 
47r e 



fJ-N + k B Ta 



A? 



^ fc2 E^E^cg 7 c(^ 7 c)+0(fc 4 )- (9) 



Here, pc and pjy are charged and neutral species bulk 
densities, pn indicates the dipole moment of a neutral 
molecule, and a at is the molecular polarizability. The 
final term sums over the the average of the square of 
given bond lengths l ai c in a charged molecule. For larger 
CHARMM- or AMBER-like molecular models, a generaliza- 
tion of this approach leading to related moment-like con- 
ditions is possible. 

Using this small-moment expansion in Eq. ([7]) and not- 
ing that the integrals of Gaussians involved can be ana- 
lytically evaluated, we find 



El 
v 



7= P cq c 



2 k R T e - 1 



+ 3 J /jp J2 PC Y1 IocCQjC (Cyc) ■ ( 10 ) 
C Q,7 

In particular, for bulk SPC/E water, which is neutral 
and nonpolarizable, we have 
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where p w is the bulk density, e w is the dielectric constant, 
and p w is the dipole moment of SPC/E water. In con- 
trast to the expression developed for ionic solution o 18 ' 19 , 
the energy correction now incorporates a significant, non- 
trivial contribution from the dipole moment. 

In fact U q /N is negative and may be bounded from 
above as Uf/N < -118.3 ^j-A 3 /cr 3 for T = 300 K, by 
assuming e — > oo and using p w determined from the rigid 
geometry of the SPC/E water molecule. In obtaining this 
numerical expression, recalling that Eq. (|10p was derived 
using cgs units is crucial. Since water has a large dielec- 
tric constant and the dipole moment contribution is large 
in magnitude, this is actually a very tight upper bound. 
If instead we use the experimental value of e w = 78, we 
find Ul/N = -118.4 ^-A 3 /cr 3 for T = 300 K with vari- 
ation lying within error bars of the simulation calculation 
of Uq. In Fig. Q]we used the infinite dielectric constant 
in our calculation of the energy correction. 

As seen in Fig. [1] the inclusion of this correction brings 
all of the energies from Gaussian-truncated simulations 
much closer to the Ewald energy, shown as a horizontal 
line. All energies now lie well within 1% deviation from 
the Ewald energy, some with substantially less error than 
that, whereas only the three larger a- values without Uf 
would lie within the less stringent 5% deviation suggested 
as sufficient in Ref. [ill 

For solutions of charged particles, previous researchers 
obtained similar correction terms for energies, though 
their physical basis was less transparent)^ 3 ^. The cor- 
rections by Hummer and coworkers relied on an analogy 
with the self-interaction in Ewald summations. The cor- 
rections by Wolf and coworkers drew upon the known lim- 
iting behaviors for charged fluids based on the Stillinger- 
Lovett moment conditions. But the necessary exten- 
sion to mixed charged and polar molecular systems was 
not appreciated. Combining thermodynamic perturba- 
tion theory with an examination of moment conditions 
for molecules, as in this paper, clarifies the general prin- 
ciples involved, and immediately leads to substantially 
improved energetics with a simple analytical energy cor- 
rection. 



V. ANALYTICAL PRESSURE CORRECTION 
FOR NVT AND NPT SIMULATIONS 

Deriving a similar correction for the pressure may seem 
more problematic, since the pressure cannot be exactly 
expressed using only site-site distribution functions 31 , 
and to our knowledge no analytic pressure corrections 
have ever been suggested. However, LMF theory pro- 
vides a general perspective that allows us to arrive at 
simple pressure corrections as well. To that end, we ex- 
press the pressure thermodynamically as 



P = T 



dS_ 
dV 



dU 
dV 



(12) 



T,{N M } ^ UV ' T,{N M } 

Since our Gaussian-truncated system with purely 
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FIG. 3: Plot of pressure without (Po shown in red crosses) and 
with the long-range correction (Po + Pi shown in blue circles) 
for the full range of truncation scales a studied. Error bars 
are shown for the data points, and error bars on the Ewald 
pressure are indicated by the thin dashed lines above and 
below the thick horizontal line at 0.044 katm. 
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FIG. 4: Plots of the volume per particle determined from 
NPT simulation using vo(r). Error bars are as in Fig. [3] 
Applying the P* correction as an external pressure brings 
particle volumes in the mimic system into better agreement 
with the Ewald result. 



short-ranged interactions reasonably captures the local 
order and structural variations expected to dominate the 
entropy, we expect that S ~ So to a very good approx- 
imation, and we use this in the first term on the right 
in Eq. (fl~2|). Corrections from the long-ranged part of 
the Coulomb interactions to the pressure Po obtained di- 
rectly from the truncated model simulation then arise 
from the second term and are simply related to a partial 
derivative of Uf with respect to volume. 

Using Eq. (fT0|) to express U\ in terms of {Nm}, V, and 
T, we find for rigid molecules that only the contribution 
due to dielectric shielding depends on volume. Therefore, 
regardless of the site composition of the rigid species, we 
have 

P*--(™£\ ksT (13) 

This correction term is purely negative, just as we 
can deduce from our simulated Po for water shown in 
Fig. [3] Using the experimental dielectric constant of wa- 
ter, e w = 78, we find Pf = -3.624 katm-A 3 /cr 3 - As 
shown in Fig. [3l including brings nearly all pres- 
sures into agreement with the Ewald result. For flexi- 
ble molecules, ix 2 N , (l 2 ), or could have a volume de- 
pendence leading to a contribution to P'. However, for 
dense and relatively incompressible systems, such a con- 
tribution is likely quite small. 

This analytical pressure correction also proves useful 
for NPT simulations of the Gaussian-truncated water. 
Shown in Fig.[4]is the volume per particle calculated dur- 
ing NPT simulations carried out at 300 K and 1 atm. The 
spherically-truncated water simulations with P = 1 atm 
have a higher volume per particle than the Ewald results 
as found in Ref. 0- However when Gaussian-truncated 
water is simulated with a corrected external pressure ad- 
justed to be P cx t = P) = 1 atm — P^T, c), the average 



volume per particle agrees quite well with Ewald results 
for all a but the smallest of 3.0 A. The latter discrepency 
simply indicates that the second order fc-space expansion 
for p qq is insufficient for the smallest a used. 



VI. CONCLUDING REMARKS 

In general, as has been well establishe d 18 ' 19 , despite 
the highly accurate local structures obtained when using 
reasonable spherical truncations, the impact of the long- 
ranged forces on thermodynamics cannot be neglected. 
We have shown here that high accuracy is possible for en- 
ergy, pressure, and density in spherically-truncated sim- 
ulations of bulk molecular fluids solely by using simple, 
analytical corrections. 

Thermodynamic corrections for nonuniform systems 
treated via LMF theory will be less straightforward. For 
example, for many of the slab systems we have simu- 
lated, a self-consistent Vr is crucial for the structu re 15 ' 16 , 
and the full LMF theory should be used for the thermo- 
dynamics as well. Corrections to thermodynamics from 
a strong-coupling simulation could perhaps be found in 
some cases based on the Carnie-Chan sum rules^, a 
nonuniform analog of the Stillinger-Lovett moment con- 
ditions, but further theoretical development is neces- 
sary. However the simple analytical corrections presented 
herein should be immediately useful in correcting the 
thermodynamics of many bulk systems of interest. 

This work was supported by NSF grants CHE0517818 
and CHE0848574. JMR acknowledges the support of the 
University of Maryland Chemical Physics fellowship. 
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APPENDIX A: DERIVATION OF ZEROTH AND 
SECOND MOMENT CONDITIONS FOR A 
MIXTURE OF NEUTRAL AND CHARGED 
SITE-SITE MOLECULES 

In this Appendix we define the two-point intermolec- 
ular charge-density function p qq (r,r') used to determine 
the total Coulomb energy and connect its behavior at 
small wave vectors to that of the fundamental charge- 
charge linear response function used in the theory of the 
dielectric constant, as discussed in Appendix B. These re- 
sults are used in the main body to derive corrections to 
the thermodynamics of uniform site-site molecular mod- 
els simulated with spherically truncated Coulomb inter- 
actions. 

The total Coulomb energy obtained during simula- 
tion of a small site-site molecular species without any 
intramolecular charge-charge interactions is 



N M N m 



U q 



M M< i=l j=l 
X 



a2i E> I i 

a=l 7=1 r iM 



q a Mq~/M' 



1 jM> 



(Al) 



In this notation, the angular brackets indicate a nor- 
malized ensemble average, M and M' indicate a given 
molecular species, i and j indicate a given molecule of 
a given species, and a and 7 represent the intramolecu- 
lar site o 23 ^ 28 . The Kronecker deltas are necessary to ex- 
clude any charge-charge interactions between intramolec- 
ular sites within a given molecule. This energy U q can 
be more compactly represented as 



dr / dr 



,p qq (r,r') 



(A2) 



where p qq is a two-point intermolecular charge-density 
function that explicity excludes any purely intramolecu- 
lar charge correlations, as implied by Eq. (|Alj) and de- 
tailed below. 

Comparing Eqs. (|A1[) and (|A2[) . we see the composite 
function p qq (r,r') is a charge- weighted linear combina- 
tion of all intermolecular, two-point site-site distribution 
function a 23 ' 28 : 

p qg (r,r') = q a Mq 1 M'PaM 1 M'iy^')- (A3) 

aM 7M' 

For our purposes here, it is more useful to relate this func- 
tion to the basic charge-charge linear response function 
used in the theory of the dielectric constant. 

For solutions of primitive model ions, Stillinger and 
Lovett showed that charge neutrality and screening place 
specific requirements on the behavior of the charge den- 
sity in fc-space at small wave vectors 2 ^. More generally, 
for a fluid composed of charged and polar molecules, the 
dielectric screening behavior of the molecules places re- 
strictions on the decay of the two-point charge density 



p qq . Based on this observation, we are able to harness a 
theoretical development of Chandler 2 - that expresses the 
dielectric constant e of polar molecules in terms of an ex- 
act sum of charge-density-weighted site pair correlation 
functions. We generalize the derivation to include both 
charged and neutral site-site molecules and we take the 
dielectric constant as a given. From this vantage point, 
we may instead use these relations to place requirements 
on the decay of the two-point charge density p qq . 

We first define the instantaneous single-point total 
charge-density p q {r, R), a function of both a given ex- 
ternal position r and the set R = {Rjjyf} = j of 

positions of all mobile charged sites in a given configura- 
tion, as 



N M nm 



(A4) 



M 1=1 a=l 



With such a definition, the ensemble-averaged charge 
density profile p q (r) is 



p q (v) = (p q (r,K)) 



(A5) 




In the case of a uniform system, p q (r) = 0. 

Comparing Eqs. (|ATj) and (fA^jl and using Eq. l(A4j) . 
we may also express the two point charge function 
p qq (\r — r'|) for a uniform system as 

p qq (\r-r'\) = (p q (r,R)p q (r / ,R)) 



(A6) 

We have used the equivalence of all molecules of type M 
in the last term. This term removes purely intramolec- 
ular charge-density correlations; we shall determine the 
small-A: contributions from this term based on well-known 
molecular properties using the approach of Chandler 23 
later in this appendix. 

The first term, in contrast, is exactly the charge-charge 
linear response function for a uniform neutral system: 

(p«(r,R>V,R)) = (<^(r,R)W,R)) 

= X qq (\r-r'\). (A7) 

Here <5p«(r,R) = p q (r,K) - (p q (r,K)). Physically X qq 
describes the coupling between charge-density fluctua- 
tions at positions r and r'. As is well establishe d 20 ' 23 , 
such a function is intimately related to the dielectric be- 
havior of the fluid at long distances, and furthermore, 
may be easily analyzed based on basic electrostatics and 
standard definitions of the functional derivative. This 
relationship is discussed in Appendix [B] 

Our goal is to write a small-fc expansion of the two- 
point intermolecular charge density, 



p qq (k) w p {0)qq + k 2 p {2)qq + 0(k 4 ), 



(A8) 
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where from Eqs. (|A6[) and (|A7|) 
^(|r-r'|) = X «(|r-r'|) 

/ nil UK 

~ ( E ^ E E 9"A^7A^( r - riM^O"' - r l] 

\ M a=ly=l 

(A9) 

As shown in Appendix [B] the charge-charge linear re- 
sponse function x qq may be expanded as 

X qq (k) = 0+^^l-^jk 2 + O(k i ). (A10) 

Now we must remove the intramolecular contributions 
as described by the last term in Eq. (|A9[) . Defining the 
conditional singlet intramolecular site density functions 
£>a|7Af(r|r') for a ^ 7 as 

P 7 M(r'K| 7 M(r|r') = (n m 5(t - r& )5(r> - r<£ )) , 

(All) 

and applying consequences of uniformity, Eq. (|A9[) can 
be written as 



^(|r-r'|)=^(|r-r'|) 

- ^2 Pm gQAfg 7 M^Q 7 M(|r - r' 

M a, 7 



(A12) 



where 



w Q7 M(|r - r'|) = 5 ai 5{v - r') + Q a \^ M {\* - r 'D- ( Al 3) 

For neutral molecules, Chandler demonstrated that the 
small- fc components of u) a ~iM{k) are related to simple 
properties of the molecule. For both charged and un- 
charged molecules, the zeroth moment of ui ai M is simply 

= S Q1 + (1 - S ai ) = 1. (A14) 
Using this exact expression in Eq. (|A12|) yields 



p(0)qq _ ^(0)qq 



E 



E P M E 1 aM( li M 

M 11,7 M 

(A15) 

an expression encompassing the standard zeroth moment 
condition for ions^l and the zeroth moment for neutral 
molecular species^. 

The expression for cD^ determined by Chandler— may 
be written most generally as 



-(2) 



"(2) 



q^7 



= ~K E q a Mq 7 M {l 2 aiM ) > ( Al6 ) 



q^7 



where l ai M is the bondlength between sites a and 7 for a 
molecule of species M. As shown in Ref.l23l for a neutral 
molecule indicated by N below, the final summation in 
the above equation is simply related to the molecular 
dipole moment fijy and the molecular polarizability a at 
as 



/ 2 



■,(2) 



1 



N 2-^i QaNq-fNU^N - 
a#7 



keTa 



JY- 



(A17) 



This relationship does not hold for a charged molecule 
since the dipole moment then depends on the choice of 
coordinate system. 

Distinguishing charged species (C) and neutral species 
(N) where {M} = {N} U {C}, and without substituting 
we find 



for 4 2) and 



P 



l2l „„ ^e-l^^^ 



An e 



- (2) 



M ct^7 



k R Te-l ^ „( 2 ) r-s .(2) ,. 10 \ 
2_^p N uj y N '-^ Pc uJc '■ (A18) 



47r e 



A' 



c 



Thus, we may write a general expression for p qq in k- 
space. Utilizing the expressions for tiff and u>q , we 
have 



P qq {k) 



E 



pcq c 



,k B Te - 1 
47r e 



k 2 ^2p N j^Mw +k B Ta A r| 

fc2 jE^E^7c(^c)+o(fc 4 ) 



(A19) 



Unlike the simple expansion of x qq (k) in Eq. (|A10|1 . we 
see that the small-fc behavior of p qq {k) depends on sev- 
eral simple properties of the solution as a whole, like the 
dipole moment and polarizability of individual neutral 
molecules, and the net molecular charge and the aver- 
age square bond lengths of charged molecules, as well 
as the dielectric constant. Thus by knowing simple sin- 
gle molecule properties and the long wavelength dielec- 
tric constant, we know how intermolecular charge-charge 
correlations decay in solution. This is the essential idea 
used to develop energy and pressure corrections for simu- 
lations of bulk liquids using molecular models with trun- 
cated Coulomb interactions. 

A related expression may be developed for larger 
molecular species with intramolecular charge-charge in- 
teractions given by CHARMM- or AMBER-like molecular 
models by modifying the total Coulomb energy to solely 
exclude the charge-charge interactions of sites i and j 
within two bonded connections of one another, using 
a "bonding function" BM{i,i) that acts similar to the 
product of Kronecker deltas in Eq. (|A1[) . In such cases, 
the expansion of x qq remains the same but the remaining 
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contribution to p qq may no longer be simply represented 
using only whole-molecule properties such as the dipole 
moment and polarizability. 



The total electrostatic potential at position r in the 
nonuniform fluid is then given by the sum of the external 
potential and the induced polarization potential: 



APPENDIX B: EXACT MOMENT CONDITIONS 
ON CHARGE-CHARGE LINEAR RESPONSE 

We present the following moment conditions for x qq 
distinct from the molecular-specific analysis found in Ap- 
pendix [A] because the results are more general than the 
specific site-site molecules chosen. The analysis of the 
behavior of x qq at small k is similar to that found in Ref- 
erences and and the final results are identical. Our 
derivation is simpler because we focus directly on the to- 
tal charge density; this also allows us to derive both the 
Stillinger-Lovett moment conditions for charged systems 
and the formula for the dielectric constant of a polar mix- 
ture from the same footing. 

The electrostatic potential at r induced by a fixed ex- 
ternal charge distribution p* xt (r') (e.g., a test charge Q 
placed at the origin, as considered by Chandler—) is given 
by 



r - r 



-dr', 



(Bl) 



and the associated electrostatic energy for a particu- 
lar microscopic configuration characterized by the set of 
molecular positions R is then 



^xt(R) 



p 9 (r,R)V oxt (r)dr. 



(B2) 



Here p 9 (r,R) is the total configurational charge density, 
defined in the particular case of a mixture of site-site 
molecules by Eq. (|A4[) . This energy contribution will 
appear in the nonuniform system's Hamiltonian when 
V e xt( r ) is nonzero. 

As such, we know from standard definitions of func- 
tional differentiation of free energie s 20 i 32 that 



S[-0A] 
5[-/?Voxt(r)] 



= <p*(r,R)) v = /$(r), (B3) 



where (3 = (keT) -1 and the subscript V indicates that 
the ensemble average is taken in the presence of an ex- 
ternal potential. Similarly we have 



H( r ) 



6 [-(3 A] 



S [-/3V cxt (r')] S [-/3Vext(r)] 5 [-/3V cxt (r')] 

= X??(r,r'). (B4) 



V to t(r) = V cxt (r)+V po i(r) 

= V cxt (r) + Idv'^L 



(B5) 



To get a formula for the dielectric constant we expand 
about the uniform neutral system and evaluate V po i to 
linear order in V ex t using Eq. (|B4| . This gives 



V po i(r 



dr' 



r r 

dr' 



r — r 



dv"P X 9q (\r'-r"\)V ext (r"). 

(B6) 



Here x qq is the linear response function in the uniform 
fluid as in Eq. (|A7p . Taking the Fourier transform of the 
final equation, we find 



V tot (k) = V oxt (k) - ^p x qq (k)V CKt (k). (B7) 



Thus to linear order we have 



Vext(k) k 



(B8) 



Phenomcnologically, we know that in the limit of k — > 
0, this ratio of the total electrostatic potential to the 
externally- imposed potential is exactly 1/e. Therefore, 
we find for our molecular mixture the general result 



1 



(B9) 



Based on the limit in Eq. (|B9[) , and expanding X qq f° r 
small k as x {0)qq + x^^k 2 , we have 



x (0)qq = 

4TT[3x (2)qq = 1- -. 



(BIO) 



Any mixture with mobile ions acts as a conductor with 
e = oo in Eq. (|B10[) . independent of the nature of the 
neutral components, and these results then reduce to the 
Stillinger-Lovett moment conditions 2 ^. 



* Present address: Materials Science Division, Lawrence 

Berkeley National Laboratory, Berkeley, CA 94720 
t Electronic address: jdw@ipst.umd.edu 
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